It is not possible, or necessary for the purpose of this investigation, to load and constrain the cranium to replicate muscle action, so simplified boundary conditions were applied and physically loaded onto the skull. The cranium was constrained by resting it on the two mastoid processes and the meso-lingual cusp of the first left molar (Fig. 1). A total compressive vertical force of 750 N was applied to the frontal bone, 15 mm anterior to bregma in increments of 75 N using a universal material testing machine (Lloyd's EZ50, Ametek-Lloyd Instruments Inc., UK). The mastoids acted as a fulcrum at which the skull rotated upon application of the downward compressive force to the frontal bone. The tooth acted as the biting point, transferring the load onto a 5-kN load cell (Omega DP25B-S; Omega Engineering Inc., Stamford, CT, USA) previously calibrated and repeatedly checked and recalibrated between experimental runs, by applying known compressive loads with the Lloyd's testing machine described above. The cell was used to measure the reaction force at the tooth during loading and was used to check it scaled linearly with the compressive force applied by the testing machine.

A Q100 DSPI measuring system (DANTEC Dynamics GmbH, Ulm, Germany), which provides a maximum field of view of 25 × 33 mm, was used to measure the strains experienced by the infra-orbital plate during mechanical loading. That region was covered with a thin layer of white spray (DIFFU-THERM developer BAB-BCB; Technische Chemie KG, Herten, Germany) to prevent light reflection. The DSPI sensor was then attached to a tripod that was glued to the infra-orbital plate using an acrylic-based adhesive (× 60; HBM Inc., Darmstadt, Germany). The software Istra Q1002.7 (DANTEC Dynamics GmbH, Ulm, Germany) was used to characterize the surface of the region of interest before loading (Fig. 1B) and create phase maps prior to and during loading (Fig. 1A) that were used to estimate maximum (ɛ1) and minimum (ɛ3) principal strain magnitudes and directions. The software then converts the strains into 2D (Fig. 1D) strain contour plots. A 3D contour plot surface of the measured maximum and minimum principal strains is also created that may be superimposed onto the surface of the virtual model (Fig. 1E). In each experimental round, the cranium was loaded with a total of 750 N in 10 increments of 75 N, as previously described. In between rounds, the load was removed and the cranium was checked to ensure it did not move and that the load was completely removed. From a total of eight experimental rounds, four were inconsistent in strain magnitudes and/or directions, attributable to experimental errors (in reading consistent strain magnitudes, possibly due to undetected interferences/slippage; and, eventually, bone failure due to the high number of experimental rounds, after which no further experimental loading was performed). Thus, these four rounds were discarded and only the consistent rounds are presented and used as reference for the FEA simulations.

The FE model was based on a CT scan of the cadaveric head, performed before dissection and mechanical testing, using a Siemens 16-channel multidetector CT scanner equipped with a STRATON tube (Siemens Somatom Sensation 16; Siemens Healthcare, Erlangen, Germany) at 120 kV and 320 mA with an H60s edge-enhancing kernel, with an original voxel size of 0.48 × 0.48 × 0.7 mm that was later resampled to 0.349 × 0.349 × 0.349 mm. Three materials were orgininally segmented by Toro-Ibacache et al. (2016), cortical bone, trabecular bone and teeth. Further segmentation was performed by thinning cortical bone in regions that presented cancellous bone. The cortical bone segmented from the internal surface of the cortical shell was subsequently allocated to cancellous bone, which was segmented as a bulk material. This additional segmentation was performed to make the model less stiff because preliminary FE simulations yielded strains that were much lower than, but proportional to, those from the experiment. Moreover, because the resolution of the scan is relatively low, partial volume averaging artificially thickens cortical bone, thus requiring further segmentation in the cortical-cancellous bone boundaries. The refined model was later exported as a BMP image stack and directly converted into a voxel based finite element mesh using the bespoke software vox2vec. Pre- and post-processing used the custom FEA program VoxFE (Fagan et al., 2007 and Liu et al., 2012). Isotropic material properties were allocated to all materials. Young's modulus of cortical bone was calculated by nanoindentation (cortical bone 16.3 ± 3.7 GPa – 21.9 ± 2.7 GPa) as previously reported by Toro-Ibacache et al. (2016) and within the range of previously reported values (Dechow et al., 1993 and Schwartz-Dabney and Dechow, 2003). Trabecular bone was allocated an elastic modulus of 56 MPa, according to Misch et al. (1999). Teeth were modeled as single strucures with a modulus of 50 GPa, within the published values for enamel (for a review, see He et al., 2006 and He and Swain, 2008). All materials were allocated a Poisson's ratio of 0.3. Model segmentation and material properties of this reference model (model 3) were then simplified by creating a two materials (cortical bone and teeth; model 2) and a one material model (cortical bone; model 1). In model 2, trabecular bone was combined with cortical bone. In model 1, trabecular bone and teeth were integrated into and represented as material with the properties of cortical bone.

The model was virtually positioned similarly to the real cranium so that the force vectors of the compressive applied loads matched between the real and virtual structures. Constraining of the model replicated the constraints applied to the real cranium. Thus, constraints were applied to both mastoid processes and to the meso-lingual cusp of the first molar. The mastoid processes were contrained at 32 nodes in all directions (left-right, anteroposterior and superiorinferior). The first molar was constrained at 21 nodes in the superiorinferior direction.

Comparison between real and predicted strains is only meaningful if the same loading regimen is applied because differences in force magnitude and/or direction would potentially result in different strain distributions, magnitudes and directions. Thus, it was necessary to verify the force vectors had the same orientation in the real and virtual experiments. To that end reaction, forces registered by the force transducer placed under the tooth were compared to the reaction forces calculated at the tooth of the FEM. Congruence of results can be interpreted as indicative of equivalent orientations of applied force vectors.

After comparison of tooth reaction forces, the deformations of the loaded crania were compared. In vitro and predicted strain directions and magnitudes were qualitatively assessed by visually comparing strain contour plots of maximum and minimum principal strains and strain vector directions. Quantitative assessment of differences between measured and predicted strain magnitudes relied on plotting magnitudes extracted along two equivalent lines traced in the DSPI and FEA outputs (Fig. 1 and Fig. 2). This approach was chosen because the DSPI and FEA outputs differ in number of dimensions (2D in the DSPI and 3D in FEA) and in resolution, thus precluding direct (point to point) comparison between results. Line 1 was placed using vertices of the field of view of the DSPI to maximize reproducibility and crossing the whole FOV to maximize the number points at which strains were sampled. It extends from the alveolus above the lateral incisor to approximately above zygomaxillare. Line 2 was placed so it would sample a different region and extends from close to alare (one of the vertices of the FOV) to approximately the maxillary root of the zygoma. The lines traced in the DSPI were reproduced in the virtual cranium by manually best fitting the 3D strain contour plot created by the Istra Q 1002.7 program onto the virtual cranium and using it as a reference to trace the equivalent lines. The strain magnitudes extracted along each of the two lines from the experimental rounds were arithemeticaly averaged and ± two standard deviations were computed. This provided a range against which FEA results were compared.

The impacts of varying segmentations (material properties) of models were first assessed by visual inspection of the strain contour plots. Quantitative assessment of sensitivity is based on maximum (ɛ1; predominantly tensile) and minimum (ɛ3; predominantly compressive) principal strain magnitudes along the above described lines (Fig. 2A) and at thirty two points over the face (Fig. 2B). These magnitudes are compared to quantify the impacts of simplifying material properties on principal strains over the face.

To assess the sensitivity of global modes of deformation to segmentation simplifications, a principal component analysis (PCA) of changes in model sizes and shapes due to loading was performed using 51 landmarks placed over the whole cranium (Fig. 2C). An additional ‘incisor bite’ was simulated to allow the magnitudes of ‘error’ due to varying modelling parameters to be calibrated in relation to the magnitude of differences resulting from application of a different bite point. The ‘incisor bite’ used the reference model 3 and differs only in that the constraint at the molar was replaced by a similar constraint at the left central incisor.

Reaction forces measured during the consistent experimental rounds at the tooth were averaged to be used as reference for the FEA models. The mean reaction force is 424 N and that of the FEA is 425.61 N. Thus, the difference is negligible (0.38%) and the loading regimens considered equivalent.

When comparing the strain contour plot of the DSPI with the FE model, it is apparent that the virtual models are stiffer, despite presenting comparable strain patterns (Fig. 3). Moreover, there is a slight mismatch in the anatomical regions where the highest/lowest strains are experienced. This is visible, e.g., on the alveolar bone overlying the root of the left canine. In the 3D contour plot of the DSPI, the alveolar bone experiences high strain in the root and low strain in the inter-radicular region. In the FEA strain contour plot, the alveolar bone over the root of the canine shows low strain and high strains are located in the inter-radicular regions. Despite this difference, comparison between the DSPI output within the corresponding area of the FEA results shows some similarities (Fig. 3 and Fig. 4). In ɛ1, these consist of a region of strain peaks that generaly envelops the left canine and premolar, rising from the inter-radicular area between the lateral incisor/canine and the premolar/first molar and converging medially to the infra-orbital foramen. This high strain drops in the region overlying the roots of the premolar and, especially, of the canine. Above the first molar, at the root of the zygomatic process of the maxilla, there is a low strain area, more distinct in the real specimen, from which strains increase both anteriorly and posteriorly. In ɛ3, the main difference between the DSPI output and the FE model is in the lower relative strains displayed by the real specimen around the first molar and above the canine. In the remaining areas, there are similarities between the models, with a relatively highly strained region that extends from the root of the zygomatic process to the inferior lateral border of the nasal aperture. In this region, there are also foci of low strains, both in the real and virtual specimens, such as in the maxillary zygomatic process, the area next to the infra-orbital foramen and over the premolar roots (Fig. 3 and Fig. 4). The numbers in Fig. 4A indicate the anatomical points at which strains peak in the real cranium (top row) and the FE model (bottom row); they are not exactly antomically equivalent and so indicate a mismatch between the model and real specimen.

Fig. 4B shows some similarities but not consistency in detail between DSPI and FEA strain vector directions (note the vectors from FEA are automatically scaled by the software according to magnitude while those from DSPI are not). Similarities in ɛ1 are generally found in the bottom half of the strain vector and contour plots, with differences being more pronounced in the top half of each. In ɛ3, similarities are found in the left half of the plot, where supero-inferior compression (vectors oriented obliquely) predominates. Differences in vector directions are found over the root of the canine (bottom centre), both in ɛ1 and ɛ3.

In order to quantitatively compare real and predicted strain magnitudes, these were extracted along lines 1 and 2 (Fig. 2). Fig. 5 presents a plot of the extracted strain magnitudes with peak strains, labelled according to Fig. 4, where the anatomical locations of these in the model and experiment are shown. In line 1 (Fig. 5A and B), there are differences in the absolute values of ɛ1 between the predicted and the real strains. These differences are most marked throughout approximately the first third of the line. Despite these absolute differences, the general strain pattern is comparable. Thus, DSPI and FE results present the same strain peaks around the canine (points 1 and 2), despite a slight mismatch in the exact location of the peaks. When considering the second and third peaks (points 2 and 3), both DSPI and FEM outputs almost match in absolute values, with the main difference being that there is a significant drop in the experimental strains between points 2 and 3 in Fig. 5A, which is absent in the predicted strains. In ɛ3 (Fig. 5B), there are also differences in the absolute values, with predicted strains being lower than real strains. Despite these differences the overall pattern of ɛ3 (Fig. 5B) between real and predicted strains is similar, with a generally flat and regular line. The main differences between DSPI and FEA outputs are in the peak of predicted strains next to the premolar (≈ 15 mm), which is absent in the DSPI output, and the flatness of the first third of the line of the predicted strains (over the canine root), which is more irregular in the real strains.

In the scatter plot of ɛ1 from line 2 (Fig. 5C), there are two main differences between real and predicted strains. The first is the significant drop in the real strains at ≈ 7–8 mm, that anatomically corresponds to the area above the apex of the canine. In the FE model, these low strains are directly over the tooth root, not extending upwards. The second difference is the drop in the predicted strains at about the apex of the premolar (≈ 20 mm, just before point 5 of Fig. 5C). Despite these two local differences, the predicted absolute strains are mostly within the range of the real ones, with a generally similar pattern between the two. In ɛ3 (Fig. 5D), both DSPI and FE extracted strains present a remakable match. The only difference is found in the higher relative strains experienced by the FE model just next to the bite point (point 7, Fig. 5D).

The effects of simplifying the segmentation, and so reducing the numbers of different materials, are most marked near the teeth. This effect is, however, much smaller with increasing distance from the teeth. This is apparent is Fig. 5 and Fig. 6. The former samples points directly over or close to teeth and shows increasing discrepancies between the FEA models as points come to lie directly over or very close to the teeth. In Fig. 6, this effect is even clearer, with almost complete overlap between models at points throughout the face that are distant from teeth and divergence between models at points over or nearby teeth (points three, eighteen and especially ten). These findings indicate that differences in material properties allocated to teeth and bone have marked local effects on predicted strains.

Analysis of large scale deformations (Fig. 7) shows that simplification of model 3 (to models 2 and 1; two then one material) impacts on magnitude (distance from unloaded model) and a little less so on mode of deformation (divergence between vectors from unloaded to loaded models). Models 1 and 2 plot directly over each other. Thus, the effect of reducing materials from 2 to 1 is minimal but reducing from 3 to 2 impacts the mode and magnitude of deformation. Anatomically, the difference in mode of global deformation between model 3 and the two simplified models corresponds to a greater degree of torsion of the face (see inset warpings and transformation grids) towards the working side, which is seen as a relatively more superior deflection of the teeth and palate and greater supero-inferior shortening of the left orbit in model 3. These results are consistent with the graphs of facial strains (Fig. 5 and Fig. 6) that show higher strains in model 3 and greater differences at points at or near the alveolus. The impact of model simplification on mode and magnitude of deformation is small compared to the impact of varying bite point (molar to incisor).